Vortex, skyrmion and elliptical domain wall textures in the two-dimensional Hubbard 

model 
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The spin and charge texture around doped holes in the two-dimensional Hubbard model is 
calculated within an unrestricted spin rotational invariant slave-boson approach. In the first part 
we examine in detail the spin structure around two holes doped in the half-filled system where we 
have studied cluster sizes up to 10 x 10. It turns out that the most stable configuration corresponds 
to a vortex-antivortex pair which has lower energy than the Neel-type bipolaron even when one 
takes the far field contribution into account. We also obtain skyrmions as local minima of the 
energy functional but with higher total energy than the vortex solutions. Additionally we have 
investigated the stability of elliptical domain walls for commensurate hole concentrations. We find 
that (i) these phases correspond to local minima of the energy functional only in case of partially 
filled walls, (ii) elliptical domain walls are only stable in the low doping regime. 

PACS numbers: 71.27.+a,71.10.Fd,75.10.-b,75.60.Ch 



I. INTRODUCTION 

The description of inhomogeneous charge and spin 
phases in the Hubbard model it a topic of current inter- 
est, mainly because it is now generally accepted that the 
high-T c compounds, at least in the underdoped regime, 
are intrinsic electronically inhomogeneous systems. A 
powerful tool for the investigation of such inhomogeneous 
electronic states is the unrestricted Hartree-Fock (HF) 
scheme which allows for the diagonalization of reason- 
able cluster sizes. Recently an extension of this approach 
based on the Gutzwiller wave function |jj has shown to 
significantly improve the HF solutions which strongly un- 
derestimate the attraction between charge carriers. In 
this paper we extend the approach of Ref. [Q to include 
transversal spin degrees of freedoms which permits the 
description of coplanar and three-dimensional inhomoge- 
neous spin textures. 

Among coplanar spin structures homogeneous spiral 
solutions have been studied by a large variety of methods 
(see e.g. [|]||) which show that a small amount of holes 
doped into the half- filled system leads toaQ~(l,l) spi- 
ral phase which changes its direction to Q ~ (1, 0) above 
some critical concentration. Relaxing the constraint of a 
homogeneous charge distribution the formation of copla- 
nar, vortexlike phases has been investigated in Ref. Q] 
using an unrestricted HF approach. These are configu- 
rations where the antiferromagnetic (AF) spin order ro- 
tates by multiples of 2ir around the localized hole. Due 
to this twist in the magnetization their energy increases 
~ ln(L 2 ) which implies their instability in large clusters. 

A three-dimensional spin texture which is known to 
be topologically stable is the skyrmion as a solution of 
the 0(3) non-linear a- model ||. It has been studied by 
several authors also for an AF background |t],|| mainly 
concentrating on small clusters. 

However, it is still controversal wether skyrmion solu- 



tions exist also on discrete lattices. Whereas unrestricted 
HF theory for the 2-d Hubbard model predicts the decay 
of skyrmions into conventional spin-polarons 0, exact 
diagonalization studies of a small cluster within the tJ- 
model Ref. H seem to support their existence even when 
one takes the contribution of the skyrmion far field into 
acount. Within our slave-boson approach we will show 
below that skyrmions textures are stable solutions of the 
discrete two-dimensional Hubbard model, when the sys- 
tem is doped with two holes away from half-filling. 

Connected with the discovery of static stripe order 
in La 2 M0 4 and Lai.^Ndo^Sro.j^CuCU a lot of 

work has been done in order to understand the differ- 
ent domain wall structures in these compounds [^2|-^5|. 
Whereas in the Ni-doped compounds one finds the stripes 
along the diagonal with one hole per Ni-site the charge 
and spin order in the Nd-doped cuprates is along the 
copper-oxygen bond direction with one hole per every 
second copper site only. Since half-filled horizontal walls 
are at odds with HF calculations the inclusion of correla- 
tions turns out to be a necessary ingredient for the study 
of domain walls in these systems |l4],[l5). In addition 
long-range Coulomb interactions can play an important 
role in stabilizing half- filled vertical stripes jl4j . 

Based on a Landau free-energy analysis of coupled 
charge and spin-density-wave order parameters it has 
been shown in Ref. fl6|| that within some region of param- 
eter domain walls may have a spiral component. Whereas 
this type of ordering is not observed in the Ni-oxides, 
some spiral contribution cannot be rigorously excluded to 
be present in the Nd-doped compounds fll}| . According 
to our analysis presented below, elliptical stripes are not 
stable for concentrations around 1/8 but may be formed 
in the very low doping limit. 

The rest of the paper is organized as follows: In Sec. 

II we give a detailed description of the formalism, in Sec. 

III we present the results for vortex, skyrmion and ellip- 
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tical domain wall solutions respectively, and in Sec. IV 
we summarize our conclusions. 



II. MODEL AND FORMALISM 

We consider the two-dimensional Hubbard model on a 
square lattice, with hopping restricted to nearest neigh- 
bors (indicated by the bracket < i,j >) 



H 



-* E < 

<ij>,a 



(1) 



where c^l destroys (creates) an electron with spin a at 

site i, and n i a = cj a Ci,cr- U is the on-site Hubbard repul- 
sion and t the transfer parameter. For the calculations 
in Sec. Ill we take t=l. In the following we use a spin- 
rotation-invariant form of the slave-boson represen- 
tation introduced by Kotliar and Ruckenstein in Ref . jfj) . 
The subsidiary boson fields , d^ stand for the anni- 
hilation (creation) of empty and doubly occupied sites, 
respectively, whereas the matrix 



Pi 



Pi,t ^(Pi,x ~ iPi,y) 

■J$(Pi,x +iPi,y) Pi, J. 



(2) 



represents the case of a singly occupied site. Since we 
consider the mean-field limit all boson operators will be 
approximated as numbers. Besides the completeness con- 
dition 



e? + fr(p? M p i)M ) + d\ = 1 



(3) 



the boson fields are constrained by the following relations 



tr(vpfpi) + 25 M , d? = V ct ( T/t ) 



(4) 



where 7u = i a 3 are the Pauli spin matrices and r M= o = 1. 



Then, in the physical subspace defined by Eqs. (|3j]J) 
the Hamiltonian (fil) takes the form 



H 



where 



-t 



E 

<ij> ,aa\a?. 



z ; = L i (e i p i + pidi)Ri 

L i= [(l-< 2 )l-pr Pi ]- 1/2 

B 4 =[(l-^)l-pfp,]- l/a 



(5) 



(6) 
(7) 
(8) 



The matrices Li and Rj guarantee the correct behavior 
in the limit U — ► within the mean-field approximation 
and pi = TpjT 1-1 is the time- reversal transformed of pi. 

The matrix elements of Zi can be calculated by trans- 
forming to a diagonal representation for the pi (see Ap- 
pendix) . The resulting effective one-particle Hamiltonian 



describes the dynamics of particles which upon hopping 
between sites are subjected to a modulation of their spin 
amplitude and spin direction, respectively. It can be di- 
agonalized by the transformation 



(9) 



where the orthogonality of the transformation requires 
5>k(*)*<,"(9) = ***- (10) 



Given a system with N e i particles we finally obtain for 
the total energy 



Etot — —t 

<ij> 

+ ^E d: 



<ij>,aa 1 a 2 
2 



k=l 



(11) 



which has to be evaluated within the constraints (||4|, [Io|) . 
This is achieved by adding these constraints quadrati- 
cally to Eq. ( pi] ) following the procedure already ap- 
plied in the Gutzwiller limit Jlj. The resulting energy 
functional then has to be minimized with respect to the 
fcrmionic and bosonic fields which is most conveniently 
done by using a standard conjugate gradient algorithm 
since the gradients of the energy functional can be calcu- 
lated analytically. In order not to end up in pathological 
side minimas we have generally started the minimization 
from an HF Ansatz for the amplitudes <&i iCr (fc). 



III. RESULTS 

Since in a previous publication the unrestricted 
slave-boson approximation has been already applied to 
the description of collinear spin structures, we will re- 
strict here to textures with two- and three dimensional 
spin ordering. In this section we discuss the spin struc- 
ture of vortex, skyrmion and elliptical domain wall tex- 
tures which turn out to be stable energy minima within 
our slave-boson approach. Obviously on finite lattices 
one has to use open boundary conditions in order to de- 
scribe higher dimensional spin structures and the cluster 
sizes we are considering in the following are ranging from 
6 x 6 up to 10 x 10. 

The incorporation of transversal spin degrees of free- 
dom allows for the definition of spin currents Vj = — dtS. 
The flow direction of these currents is along the bonds 
of the lattice, however, they are additionally vectorial in 
spin space and within the present approach we obtain for 
the i-th component of the spin current flowing between 
sites < nm >: 
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a 2 Z n,a 3 a 2 z ra,a l a i $171,04 V s ) 



S, = 5 , e lQR >[cos( 



(13) 



(7ia20'3(T4 k — 1 

(12) 

where t 1 are the Pauli matrices and the hopping factors 
z n aa i are defined in Eq. @. The i-th component of j nm 
can be thought of as measuring the spin-twist in the or- 
thogonal directions I, k ^ i so that the total current com- 
ponents j l nl which are plotted in the results, visualize the 
direction of maximal twist in the spin components Ik ^ i 
at lattice site n. 



A. Vortex States 

The structure of vortex solutions, where the magneti- 
zation rotates in a plane by some multiples of 2n around 
the localized holes, has already been studied in Ref. (J] 
within unrestricted HF theory. We also obtain vortex 



states as local minima of the energy functional (11), 
where the total energy is about 5% — 10% lower than in 
the HF approach, depending on lattice size and on-site 
repulsion U. However, for one hole away from half-filling 
vortex solutions are always higher in energy than the con- 
ventional Neel ordered spin polaron. Moreover, their to- 
tal energy increases logarithmically with the cluster size 
as a consequence of the twist between neighboring mag- 
netization vectors in agreement with Ref. H. 
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FIG. 1. Spin structure and spin currents for a vortex- 
antivortex pair on a 10 x 10 lattice. The Hubbard on-site 
repulsion is U = Wt. 



where <pi, (4>2) refer to the angles between x-axis and the 
vectors connecting vortex (antivortex) core and site Rj. 
The AF wave vector is denoted by Q. 

Concerning the stability of the vortex-antivortex pair 
we obtain for the 10 x 10 cluster and U = lOi a binding 
energy of E h = — 0.047i with respect to the Neel type 
bipolaron. Taking into account the far field energy, cal- 
culated within the xy-model for the spin structure Eq. 
( |l3| ) and an exchange constant of J = At 2 /U = OAt still 
results in a negative binding energy of E\ ot = -0.016L 

We note that the absolute value of E\ ot slightly de- 
creases for smaller cluster sizes since the boundary spins 
do not properly adjust to the solution Eq. (13). One can 
therefore safely conclude that vortex-antivortex pairs are 
also stable in the thermodynamic cluster limit. 



B. Skyrmions 

On a discrete 2-dimensional AF lattice the spin struc- 
ture of skyrmions, originally obtained as solutions of the 
0(3) non- linear a- model M, has the form M,p| 



1 - -4- 1 - 



ii + i* + A 2 



(-1) 



s, = (-i) 



%l + i% + A 2 
1 i 2 x + i\ - A 2 
2^ + z 2 +A 2 



(14) 



where A denotes the core size of the skyrmion and its 
center is located at i x = i y = 0. 

In order to enhance convergence we initialized our min- 
imization procedure with (non self-consistent) HF wave 
functions corresponding to the spin fields Eqs. (|l4|). 

Despite intensive search we could not obtain skyrmion 
states for one hole doped in the half-filled system. These 
solutions always converged towards a spin-polaron em- 
bedded into a collinear AF Neel state as already observed 
in Ref. f§. 



This logarithmic divergency can be compensated when 
two holes form a vortex-antivortex pair. According to our 
calculations the vortex cores are located at the center of 
diagonally next nearest neighbor plaquettes, thus sepa- 
rated by the plaquette where the two holes are localized. 
Fig. 1 shows the spin structure and spin current of such a 
vortex-antivortex pair on a 10 x 10 cluster and U = lOt. 
All the spins lie in the xy-plane which means that only 
the z-component of the spin current has a non zero con- 
tribution. The spin field of this texture can be described 
by 




FIG. 2. Charge- ((n»)) and Spin-distribution for a skyrmion 
texture on a 8 x 8 lattice for U = Wt. 
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FIG. 3. xy-, yz- and xz-projections of the spins together with the respective spin currents for the skyrmion shown in Fig. 2 



The situation changes when removing two particles 
from the half-filled system. Fig. 2 displays the charge- 
and spin structure in case of 62 particles on a 8 x 8 lat- 
tice for U — lOt. The two holes then are localized on a 
2x2 plaquette at the skyrmion center and in their vicinity 
spins show a remarkable deviation from the z-direction. 
This is more easily seen in Fig. 3 where we have plotted 
the xy-, xz- and yz- spin projections, respectively, to- 
gether with the corresponding spin currents. The xy-spin 
components rotate by 360° around the skyrmion center 
resulting in a circular spin current for j z . However, all 
current components strongly decay for sites far away from 
the skyrmion center indicating that the core size parame- 
ter A is small. Upon fitting our solutions to the skyrmion 
field Eq. (jTJ) we obtain A m 0.7 for U = lOi and cluster 
sizes 6 x 6, 8 x 8 and 10 x 10, respectively. This already 
indicates that the skyrmion state should survive the limit 
of large clusters. 

To assess the question of stability in more detail we 
have also calculated the total energy using skyrmionic 
boundary conditions H. These are defined through an 
exchange field JS(R, A) to which the spins at the bound- 
ary are coupled and S(R, A) has the form of Eq. (|l4|). For 
the exchange constant we take the strong-coupling value 
of the Hubbard model J = U 2 /U. The total energy 
then is evaluated as a function of A and in the result we 
substract the energy contribution of the exchange field. 
This energy has to be compared with the corresponding 
value of a collinear bipolaron where two holes are local- 
ized on neighbored sites within the Neel ordered system. 



The results are plotted in Fig. 4 again for U = Wt and 
three different cluster sizes. As can be seen all curves 
display a clear minimum at some value of A indicating 
the presence of a stable skyrmion solution with a signifi- 
cant lower energy than the collinear bipolaron. It should 
be mentioned that these minima for the corresponding 
one-hole doped systems are always at A = 0, i.e. the 
configuration of a conventional spin-polaron. From the 
results shown in Fig. 4 one further sees that a lowering 
of energy with respect to the AF bipolaron is already ob- 
tained for A = indicating that despite the imposed Neel 
boundaries the system has a skyrmion like core. This en- 
ergy shift increases with the system size since the central 
spins of large lattices can more properly adjust to the 
skyrmion state. From the fact that the A = results for 
the 8x8 and 10 x 10 coincide within the numerical error 
we conclude that also in the thermodynamic limit the 
skyrmion solution should survive. Also plotted in Fig. 
4 with filled symbols are the energy differences between 
bipolarons and skyrmions calculated with open boundary 
conditions. As already mentioned one obtains the same 
core parameter A for all three cluster sizes which agrees 
with the position of the minimum of the solid line (10 x 10 
system). Also this feature demonstrates that our largest 
cluster should already correctly describe the skyrmion 
structure of infinite clusters. Since the skyrmion on small 
lattices is very much influenced by the boundary condi- 
tions the minimum of the dotted curve (6x6 cluster) is 
shifted to a higher value of A s» 1. 
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FIG. 4. Total energy of a two-hole-skyrmion with respect 
to the energy of a Neel-type bipolaron as a function of the 
skyrmion core size parameter A. The boundary spins have 
been coupled to the skyrmion solution Eq. (^) via a mean- 
field exchange field. Solid line and circles: 10 x 10 lattice; 
Dashed line and squares: 8x8 lattice; Dotted line and dia- 
monds: 6x6 lattice. The full symbols mark the energies for 
open boundary conditions. U = Wt. 

In order to compare the stability of skyrmion states 
with the vortex- antivortex solutions, one has to take into 
account the far field contribution to the total energy. 
Considering again a 10 x 10 system and U = 10* the 
skyrmion binding energy with respect to the Neel-type 
bipolaron (cf. Fig. 4) is E b tot = -0.0265*. Including 
the far field energy (taking again J = 0.4*) one obtains 
E h tot — —0.0071* which is approximately half the value of 
the vortex-antivortex binding energy. 



solution whereas a — reduces the spin structure to a 
collinear 'classical' domain wall. 

We have used (non-self consistent) HF states corre- 
sponding to Eq. (|15|) as starting fields for our minimiza- 
tion for different values of a. Only vertical domain walls 
have been considered and periodic boundary conditions 
where applied in the x- and y-direction, respectively. In 
case of a completely filled domain wall (i.e. one hole per 
site along the wall) we only found the collinear solutions 
whereas coplanar structures become stable for half-filled 
walls. 
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C. Elliptical domain walls 



FIG. 5. Two possible spin structures for elliptical domain 
walls together with the corresponding spin currents. The on- 
site repulsion is U = 6* and periodic boundary conditions in 
x- and y-direction have been used. Shown are the results for 
a 9 x 8 lattice doped with 4 holes. 



The possibility that charge-spin coupling in correlated 
systems may induce the formation of elliptical domain 
walls was proposed by an analysis of the Landau free- 
energy functional for coupled charge- and spin density 
waves Jl6[ . These are coplanar spin structures where the 
spin components in the first harmonic are modulated as 



S x (r) = S e lQr cos(a) cos(qr) 
S y (r) = ±S e lQr sin(a) sin(qr) 



(15) 



The eccentricity of the elliptical domain wall is deter- 
mined by a and Q,q correspond to the wave vectors of 
commensurate AF and the domain wall periodicity, re- 
spectively. The case a = ir/4 describes an ideal spiral 



This is shown in Fig. 5 for a 9 x 8 lattice with 68 
particles where we plot two kinds of possible spin struc- 
tures which are local minima of the energy functional 
Eq. (^). Also shown are the respective spin currents 
(only z-components since the magnetization vector is 
completely in the xy-planc). Spin fields and currents 
display a quadrupled structure along the wall which is 
a necessary condition for stability within any mean-field 
approach pa]. Due to the stripe charge structure the 
current flows are more complex than for the elliptical so- 



lution Eq. (15) which predicts currents flowing in a single 
direction orthogonal to the wall. Instead we observe also 
currents along the wall (Fig. 5a) and a vortex-antivortex 
structure in Fig. 5b. 
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Although the elliptical stripes are local minima of the 
energy functional Eq. (|ll|) they are slightly higher in en- 
ergy (« 5%) than collinear domain walls. These have 
been shown to correspond to the ground state when one 
takes long-range Coulomb interactions into account [jl4| . 
However, we do not expect significant differences of a 
long-range contribution to collinear and elliptical stripe 
solutions. 

The structures shown in Fig. 5 correspond to systems 
with hole doping 1/18 and we could not obtain elliptical 
solutions for higher doping. However, one can also study 
the very low doping limit upon using open boundaries in 
the x-direction. Indeed in this case elliptical half-filled 
stripes become favored with repect to collinear domain 
walls for not too large values of the on-site repulsion U 

(< m). 



a) 



XV 



r 



b) 



- N N V 


\ 


T 




y y y 


AN A 


\ 




s 


S A A - 


■ S N A 




\ 


/ 


A s a - 


N N % 


\ 




/ 


A A.A ■ 


NN N \ 


r 




y y A - 


N N N 


\ 




/• 


A A y 


- N N V 




\ 


/ 


A A y 


■ \ \ N 


\ 




/• 


S A y - 


2 


4 




6 


8 1 




XV 






- N N "x 


\ 


T 


/ 


A N A - 


- N N X 


\ 


1 


/ 


A A A - 


AN V 


\l 






y a a h 


- N N N 


\ 






S A A 


-ANA 


\ 


I 


/ 


A N A -< 


N N N 


\ 


I 


/ 


A A A 


-ANA 


\ 






y A - 


NAN 


\ 






A A ^ 



- J J , \ 




/ " 




■ -, ■, - / 




\ ~ 




- , r - \ 




/ - 




■■-!.„/ 




\ ~ 




- J J . \ 




/ - 




■ "1 1 ' / 




\ ^ 




- r r - \ 




/ - 




, >- t » / 




\ - 





FIG. 6. The same as in Fig. 5 but now for open boundaries 
in x-direction. 

In Fig. 6 we show the same spin structures as in Fig. 5 
but for open boundaries in x-direction. The main differ- 
ence between these two figures concerns the angle of spin 
rotation O across the domain wall which is only « 3/4 
of the expected value of it according to Eq. ([l5]) . This is 
due to the fact that the system now can acquire a state 
which has zero total spin current, whereas for periodic 
boundaries the spiral component always requires a net 
flow in x-direction. 

However, since for a regular array of stripes the charge 
periodicity and the spin modulation have to be related 
by k charge = 2k spm @ the domain wall induced kink 



type spin rotation has to be supplemented by an addi- 
tional spiral field. Assuming an exponential relaxation 
of the spin rotation from 8 to ir by this spiral field gives 
an additional energy per site of AE w p s (n — Q) 2 X/L 
where p s = JS 2 is the spin stiffness, A defines the length 
scale of relaxation and L is the stripe separation. Conse- 
quently this additional energy can become small enough 
in the low doping regime so that elliptical stripe configu- 
rations are more stable than collinear solutions for small 
hole concentrations. 



IV. CONCLUSION 

Summarizing, we have presented the structure of spin 
textures which we obtained by applying an unrestricted 
slave-boson mean-field approximation in its spin-rotation 
invariant form to the 2D Hubbard model. This approach 
is suited for calculating the charge and spin distribution 
of electrons (holes) in inhomogcneous and strongly cor- 
related systems since it incorporates correlation effects 
beyond the standard unrestricted HF theory. These cor- 
relation effects turn out to be especially important for the 
interaction among holes in the 2-D Hubbard model where 
the HF approximation strongly underestimates their at- 
tractive potential (!]]. 

Including the effect of transversal degrees of freedom 
we have shown that within our slave-boson mean-field ap- 
proach two holes in the half-filled two-dimensional Hub- 
bard model are bound by forming a vortex-antivortex 
pair, oriented along the diagonal direction. This tex- 
ture has significantly lower energy than a conventional 
collinear bipolaron, also when the far field contribution 
is taken into account. Additionally we have found that 
skyrmion states can be stable even on discrete lattices. 
We attribute this to the inclusion of correlation effects 
within our approach, since unrestricted Hartree-Fock the- 
ory cannot account for skyrmions as self-consistent solu- 
tions (or local minima of energy). Indeed it has been 
also observed by the authors of Ref. || that within the 
tJ-model a semiclassical description cannot account for 
the occurence of skyrmions but that it is the quantum 
fluctuations which stabilize this texture. Considering the 
formation of elliptical domain walls it turned out that 
these structures only appear for half-filled walls in the 
low doping regime. Similar to the collinear stripes they 
are stabilized by a quadrupling of the period along the 
wall. This can be realized either by alternating on- wall 
spin currents (Fig. 5a) or by forming a vortex-antivortex 
structure (Fig. 5b). 
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APPENDIX: 

Our purpose is to calculate the matrix elements of the 
hopping matrices Zj by transforming to a diagonal rep- 
resentation of the matrix pi (Eq. (0)) according to 



7XPi)=Xi 1 PiXi = 



XT 
AT 



(Al) 



where the eigenvalues are given by Xf — + 



Pi,l) ± hy(Pi,1 ~ P^i) 2 + 2 (plx+Pl y )- The transforma- 
tion matrix \i reads as 



Xi 



0t = 



<4 a i 
pf Pi 

1 Pi,x ~ ^Pi,y 



(A2) 
(A3) 

(A4) 



Applied to the hopping matrix this transformation yields 

z t = XiT(zi)xi 
T(z ; ) = T(L i )[e i T(p I ) + T(p i )d t ]T(R i ) (A5) 

and since T(pi) is diagonal it is straightforward to eval- 
uate T(Li),T(Ri) and T(p;) as 



T(U) = 

T(R0 = 
T(pi) = 



Vl-rf?-(A+) 2 





Vl-^-(A-) 2 

o 



^l_ d 2_ (A - )2 


1 



A- 
A+ 



(A6) 

(A7) 
(A8) 
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